Ultracold Fermions in a Graphene-Type Optical Lattice 
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Some important features of the graphene physics can be reproduced by loading ultracold fermionic 
atoms in a two-dimensional optical lattice with honeycomb symmetry and we address here its ex- 
perimental feasibility. We analyze in great details the optical lattice generated by the coherent 
superposition of three coplanar running laser waves with respective angles 27r/3. The corresponding 
band structure displays Dirac cones located at the corners of the Brillouin zone and close to half- 
filling this system is well described by massless Dirac fermions. We characterize their properties by 
accurately deriving the nearest-neighbor hopping parameter to as a function of the optical lattice 
parameters. Our semi-classical instanton method proves in excellent agreement with an exact nu- 
merical diagonalization of the full Hamilton operator in the tight-binding regime. We conclude that 
the temperature range needed to access the Dirac fermions regime is within experimental reach. We 
also analyze imperfections in the laser configuration as they lead to optical lattice distortions which 
affect the Dirac fermions. We show that the Dirac cones do survive up to some critical intensity or 
angle mismatches which are easily controlled in actual experiments. In the tight-binding regime, we 
predict, and numerically confirm, that these critical mismatches are inversely proportional to the 
square-root of the optical potential strength. We also briefly discuss the interesting possibility of 
fine-tuning the mass of the Dirac fermions by controlling the laser phase in an optical lattice gener- 
ated by the incoherent superposition of three coplanar independent standing waves with respective 
angles 27r/3. 

PACS numbers: 03.75.Lm, 03.75.Ss, 37.10.Jk, 71.10.Fd 
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I. INTRODUCTION 



In 2004, researchers in Manchester isolated one-atom 
thick sheets of carbon atoms, with the atoms organized 
in a planar honeycomb structure Such a mate- 

rial is referred to as graphene and is of utmost im- 
portance in condensed-matter physics since by stacking 
it one gets the graphite structure, and by wrapping it 
one gets carbon nanotubes and fuUerenes [31 . Graphene 
is also of great theoretical interest because it provides 
a physical realization of two-dimensional field theories 
with quantum anomalies Q. Indeed, the effective the- 
ory that describes the low-energy electronic excitations in 
graphene is that of two-dimensional massless Weyl-Dirac 
fermions. In graphene these massless fermions propagate 
with about one 3GGth of the speed of light. Triggered by 
the Manchester discovery, an intense activity has flour- 
ished in the field, and continues to flourish, as witnessed 
by Refs. 0, i, i, 0, B H, for example. The reported 
and predicted phenomena include the Klein paradox (the 
perfect transmission of relativistic particles through high 
and wide potential barriers) 0, the anomalous quantum 
Hall effect induced by Berry phases [lH , and its cor- 
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responding modified Landau levels [12| . 

It is now well established that some condensed-matter 
phenomena can be reproduced by loading ultracold 
atoms into optical lattices [H, [l^. The great advan- 
tage is that the relevant parameters are accessible for 
accurate control (shape and strength of the light poten- 
tial, atom-atom interaction strength via Feshbach res- 
onances [l^, etc.) while spurious effects that destroy 
the quantum coherence are absent, such as the analog of 
the electron-phonon interaction. Our present objective 
is to analyze a scheme capable of reproducing in atomic 
physics the unique situation found in graphene [l6| . It 
consists of creating a two-dimensional honeycomb opti- 
cal lattice and loading it with ultracold fermions like the 
neutral Lithium-6 or Potassium-40 atoms. 

Parts of this paper recall known results. In addi- 
tion to the need of setting the stage and introducing 
the notational conventions, there is also the intention to 
bridge the solid-state community and the atomic physics 
community on the particular subject of massless Dirac 
fermions as observed in graphene sheets and its counter- 
part in atomic physics. We also present extensions of 
previous solid-state works in the atomic physics context 
and report a number of new results. 

We analyze the various experimental parameters that 
need to be controlled in order to reproduce, with cold 
atoms trapped in an optical lattice, the physics at work 
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in graphenc. After briefly introducing optical lattices, 
we first explain how to create an optical lattice with the 
honeycomb symmetry and analyze its crystallographic 
features. We then calculate the band structure in the 
tight-binding approximation and by exact diagonaliza- 
tion, thereby providing evidence for the occurrence of the 
so-called Dirac points. Next, we evaluate the nearest- 
neighbors hopping amplitude by using a semi-classical 
instanton method. For the benefit of possible experi- 
ments we give the necessary requirements for reaching 
the massless Dirac fermions regime. Finally, we examine 
how massless Dirac fermions survive lattice distortions 
that could result from intensity-unbalanced or misaligned 
laser beams. These distortions open the way to new 
physics related to the quantum Hall effect [l3|- We will 
close by briefly mentioning possible experiments to tar- 
get for noninteracting and interacting ultracold fermions 

MM- 

II. THE HONEYCOMB OPTICAL LATTICE 

A. Radiative forces and optical lattices 

A two-level atom (with angular frequency separation 
LUat and excited-state angular frequency width F) that 
interacts with a monochromatic laser fleld with complex 
amphtude £{r,t) — £J(r) e~"^^* gets polarized and ex- 
periences radiative forces due to photon absorption and 
emission cycles [1^ [2l|. When the light frequency is 
tuned far away from the atomic resonance, i.e., when 
the light detuning S ^ ujl — Uat is much larger than 
F, the fleld-induced saturation effects are negligible and 
the atom essentially keeps staying in its ground state. 
In this situation, the atom-fleld interaction is dominated 
by stimulated emission processes where the atomic dipolc 
absorbs a photon from one Fourier component of the field 
and radiates it back into the same or another one of these 
Fourier modes. In each such stimulated cycle, there is a 
momentum transfer to the atom and, as a net result, the 
atom experiences an average force in the course of time. 
This dipole force exerted by the field onto the atom in 
its ground state is conservative. It derives from the po- 
larization en ergy shift of the atomic levels (AC Stark or 
light shifts) [24] and the dipole potential V{r) is given 
by 

where I{r) = eoc|£?('r)|^/2 is the light field intensity 
(time-averaged energy current density) at the ccnter-of- 
mass position r of the atom and Is is the saturation in- 
tensity of the atom under consideration. 

For multi-level atoms, the situation is more compli- 
cated as the dipole potential now depends on the par- 
ticular atomic ground state sub-level under considera- 
tion. However, if the laser detuning S is much larger 



than the fine and hyperfine structure splittings of the 
atomic electronic transition, then all ground state atomic 
sub-levels will essentially experience the same dipole po- 
tential. This common potential turns out to be given by 
^ as well. Hence, by conveniently tailoring the space 
and time dependence of the laser field, one can produce 
a great variety of dipole potentials and thus manipulate 
the ground state atomic motion. 

Optical lattices are periodic intensity patterns of light 
obtained through the interference of several monochro- 
matic laser beams [2^. By loading ultracold atoms into 
such artificial crystals of light one obtains periodic arrays 
of atoms. Indeed, as seen from ([1]), when the light field 
is blue-detuned from the atomic resonance ((5 > 0), then 
the atoms can be trapped in the field-intensity minima 
whereas for red-tuned light {6 < 0) they can be trapped 
at the field intensity maxima. Such arrays of ultracold 
atoms trapped in optical lattices have been used in a 
wide variety of experiments. As recently evidenced by 
the observation of the Mott-Hubbard transition with de- 
generate gases [15 , they have proven to be a unique tool 
to mimic, test and go beyond phenomena observed un- 
til now in the condensed-matter realm [13, UHl- They 
also have a promising potential for the implementation 
of quantum simulators and for quantum information pro- 
cessing purposes [H, [2^ [13] • 

B. Optical lattice with honeycomb structure 

1. Field configuration and associated dipole potential 

The simplest possible optical lattice with honeycomb 
structure is generated by superposing three coplanar 
traveling plane waves that have the same angular fre- 
quency ujl = ckL, the same field strength Eq > 0, the 
same polarization and the three wave vectors fea form a 
trine: their sum vanishes and the angle between any two 
of them is 2tt/3, 

fci+fe2 + fc3 = 0, ka-kb = kl(^hab-^) (2) 

with a,b ~ 1,2,3 and Sab is the Kronecker symbol [23| . 
As is illustrated in Fig. [1] we choose the x, y-plane as the 
common plane of propagation and, to be specific, use 

k.-k,e,, (3) 

for the parameterization of the wave vectors. 

Further, we take all fields to be linearly polarized or- 
thogonal to the plane, so that the three complex field 
amplitudes are given by 

£air, t) = Eq S^- ■ ^ ^ <t>a)^-~i^Lt 

where 0^ is the phase of the ath field for < = at r = 0. 
We note that a joint shift of the reference points in time 



FIG. 1: The coplanar three-beam configuration used to gen- 
erate the honeycomb lattice. All beams have the same fre- 
quency, strength and linear polarization orthogonal to their 
common propagation plane. The honeycomb lattice under 
consideration is obtained for blue-detuned beams with respec- 
tive angles 27r/3. For these symmetric laser beams, the time- 
averaged radiation pressure — albeit small at large detuning 
— vanishes in this configuration. By reversing the propaga- 
tion direction of one of the lasers, such that fci = k2 + ks, say, 
a triangular lattice of a different geometry is formed. We will, 
however, exclusively deal with the fci fe2 + fcg = case. 

and space, 

a ^ a 

removes the phases (pa from ([4]), so that the simple choice 
01 = ^2 = '/'s = is permissible, and we adopt this con- 
vention. In an experimental implementation, one would 
need to stabilize the phase differences (f>a — 4>b to prevent 
a rapid jitter of the lattice that could perturb the atoms 
trapped in the potential minima. 

The dipole potential ^ generated by the electric field 
E = Sa is of the form 

V{r) = Vo\l{rf = Vov{r) with = , (6) 

where Iq is the intensity associated with the field strength 
Eo. The total dimensionlcss field amplitude f{r) and the 
dimcnsionlcss optical potential v{r) arc given by 

/(r) = 1 + exp(-i6i • r) + exp(i62 • r) (7) 

and 

v{r) = 3 + 2cos(6i-r) + 2cos(62-r) + 2cos((bi + 62) • r) , 

(8) 

where 61 = fca — fci and 62 = fei ^ ^2 feature the recip- 
rocal primitive vectors. For the parameterization ([3]), we 
have 

b2 } 2 ^ ' 

with K = 1 6a I = \/3fcL. One may further notice that the 
periodic patterns associated to each of the cosine terms 
in ^ have the same spatial period of (27r/fci)/\/3, about 
58% of the laser wavelength Al = 27r/fci. 



FIG. 2: The triangular reciprocal lattice B associated with 
the triangular Bravais lattice of Fig. [l] It is spanned by the 
reciprocal primitive vectors 61 and of ((9|, and is also a 
triangular lattice (as indicated by the full dots). The shaded 
region identifies the first Brillouin zone which is here a reg- 
ular hexagon. Its center is conventionally named F in the 
solid-state literature. Opposite edges are in fact identical as 
they only differ by a translation in the reciprocal lattice. This 
feature is emphasized by drawing the identical edges with the 
same (solid, dashed or dash-dotted) line. For the same rea- 
son, the three corners Ka {a = 1, 2, 3) are to be identified with 
each other, and likewise the three corners K'^ are really only 
one point in fi. Thus only two of the six corners, collectively 
labeled as K and K' and known as the Dirac points, are dif- 
ferent. Also shown are the wave vectors of the three coplanar 
plane waves (dashed arrows). 

Linear combinations of the Brillouin vectors with inte- 
ger coefficients define the reciprocal lattice B, a regular 
pattern in fc-space, 

B = {nibi + 71262 I ni, 712 = 0, ±1, ±2, . . . } . (10) 

The reciprocal lattice is central to all studies of the dy- 
namics of particles that move under the influence of the 
given periodic potential [28t . 

In particular, one domain in reciprocal space of utmost 
importance is the first Brillouin zone Q, defined as the 
so-called primitive Wigner-Seitz cell [11] of B, see Fig. [5] 
It is a regular hexagon but with the subtle feature that 
opposite edges are to be identified with each other since 
they can be related by a displacement vector in B. By 
the same token the three corners Ka (respectively K'^) 
have to be identified with one another and we collectively 
denote them by K (respectively K'). These two different 
corners K and K' are known in the graphene literature 
as the Dirac points for a reason that will become clear in 
the next section. Upon denoting K = Ki and K' = K'^, 
their positions in f2 are given by the wave vector of the 
lasers that generate the optical honeycomb potential, 

K ^ -K' ^ ^{b2 - bi) = (11) 

and K2 = fc2 = -fC — 62 ■ = k^ = K + bi, as well as 
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FIG. 3: The underlying Bravais lattice Z3 of a two-dimensional 
honeycomb is the two-dimensional triangular Bravais lattice 
with a two-point basis A and B. The grey-shaded area is the 
primitive cell E. The honeycomb lattice constant a is defined 
as the distance between nearest-neighbor sites. 



2. Triangular Bravais lattice 

The dimcnsionlcss potential ([S]) consists of a periodic 
two-dimensional array of maxima, minima, and saddle 
points, generated by repeated translations of a primitive 
unit tile called the basis. The underlying lattice geometry 
itself is encapsulated in the associated Bravais lattice B, 
that is 

B = {miOi + ?T!,2a2 I mi, 7712 = 0, ±1, ±2, . . . }, (12) 

such that the value of the potential is not affected by any 
displacement R G B, v{r + R) = v{r). 

The Bravais primitive vectors are constructed based 
on the relation 



b- 



(13) 



In other words, the Bravais lattice B and the Brillouin 
lattice B constitute dual spaces. Supplementing we 
have the explicit parameterization 



ai 
0.2 



= A 



V3 



(14) 



where A = \aa\ ~ 47r/(3fcL) = 2Al/3 is the common 
length of the Bravais primitive vectors. 

The Bravais lattice defined by (jl4p is a triangular one. 
We opt here for the diamond-shaped primitive cell E de- 
lineated by the two Bravais lattice vectors as a tiling 
for the optical potential ([5]); see Fig. [31 Another possi- 
ble choice would have been the hexagonal Wigner-Seitz 
cell . This cell is useful when discussing the symmetry 
group of the lattice. 

To proceed further one now needs to analyze the struc- 
ture of the optical potential ([8|) inside the primitive cell. 
In passing, we mention here that red detuned {d < 0) 
lasers give Vq < and there is only one potential mini- 
mum in each primitive cell E. Upon trapping atoms in 
these potential minima, one gets a triangular lattice that 
is not of graphene type. This situation is interesting in 
view of quantum magnetism and frustration phenomena 
but it is not the situation we want to study here. 



3. The honeycomb structure 

When the optical lattice is instead blue-detuned {6 > 
0) , Vb is positive and atoms are "weak- field seekers" . The 
potential minima coincide with the minima of the elec- 
tric field strength, and the maxima coincide as well. By 
choice of coordinate system, the maxima locate at the 
Bravais sites and the dimensionless potential ^ has its 
maximal value of v{0) = 9 at the corners O, P, Q, R of 
the diamond-shaped primitive cell E, see Fig.[4l 

Two different potential minima, given by the zeros of 
the total dimensionless field amplitude /(r), are found 
in E at 



A 

7!^ 



and 



2r, 



(15) 



respectively. From a crystallographic point of view, E 
is a primitive cell with a two-point basis. By applying 
repeated Bravais translations on E, one generates two 
different sublattices of potential minima, one made up 
of A-type sites and the other made of B-type sites, see 
Fig. [3] and Fig. [D Altogether the potential minima are 
organized in a honeycomb structure reminiscent of the 
positions of the carbon atoms in graphene sheets. 

The three displacements that move an A site to a neigh- 
boring B site — they translate the A sublattice to the B 
sublattice — are parameterized by 



1 



C2 

C3 



(ai + 02) = ae^ 



+ V3e 



-(a2-2ai) = a- ^ 
-(ai - 202) = a ^ ^ 



(16) 



where a = |cj | = A/\/3 = 47r/(3K) = 2Al/\/27 is the 
honeycomb lattice constant. It is the distance from an 
A site to a neighboring B site, or the distance from the 
center of the hexagon of minima to one of its corners. 

Halfway between two neighboring minima, the poten- 
tial has saddle points where vir) = 1. They are lo- 
cated at the center and at the middle of the edges of 
E, see Fig. ID As the saddle points on opposite sides of E 
are connected by Bravais displacements, there are there- 
fore three nonequivalent triangular sublattices of saddle 
points, and we thus count three saddle points per primi- 
tive cell. 

We also note that the potential is invariant under 120° 
rotations around the locations of the potential minima 
and maxima and, therefore, that the potential is isotropic 
in the vicinity of these points. We anticipate that the 
local harmonic oscillator potential at a minimum will be 
isotropic; see ([35)1 below. By contrast, the corresponding 
local potential at a saddle point is not isotropic. 

All these matters are illustrated in Fig. 01 where we 
clearly identify the various triangular sublattices. Cold 
fermionic atoms trapped in this optical potential would 
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FIG. 4: [Color online] Left: The honeycomb pattern composed of the triangular lattices of minima at sites A and B, of maxima 
at sites C, as well as of the saddle points between neighboring A and B sites (marked by dots). The bottom plot shows the 
potential along the x axis which is one of the . . . ABCABC. . . lines with 2: = at a C site. The saddle points S appear as local 
maxima here, with a height that is one ninth of the global maxima at sites C. Cold atoms trapped in this optical potential 
would be found at the A and B sites. Right: Equipotential lines for the optical honeycomb potential (|6}. Along the straight 
black lines that connect the saddle points, we have V{r) = Vb- The red closed circular curves fill out a hexagonal area, centered 
at the points of maximal potential; from inside out the respective values are V{r) = 8Vb, 5Vo, 2Vb, and 1.05Vb. The closed 
curves in blue and green fill out areas of the shape of equilateral triangles, their centers are the minima that constitute the A 
sublattice (blue) or the B sublattice (green); along the curves the potential has the values V{r) — 0.95Vo, O.6V0, 0.3Vb, and 
O.OSVq. One primitive diamond-shaped unit tile S spanned by ai and 02 is traced out. It contains two different minima, one 
of A-type (in blue, on the left inside) and one of B-type (in green, on the right inside). The trine of the A — > B displacement 
vectors (|16p is indicated as well. Finally, for completeness, we also trace out the Bravais Wigner-Seitz unit tile. It is a hexagon 
centered at a potential maximum and with potential minima at its corners. 



be found at the A and B sites, similar to the binding of 
electrons in graphene to the carbon ions. 

As a side remark, it may be worth mentioning that the 
saddle points affect the classical dynamics of a particle 
evolving in the honeycomb potential with a sufRciently 
large energy. Since the potential is nonseparable and an- 
gular momentum is not conserved here, the saddle-points 
could be the seed for instabilities in which case the mo- 
tion could turn out to be nonintegrable and chaotic. If 
so, this chaotic behavior should then be revealed, for ex- 
ample, in the statistical properties of the quantum spec- 
tra, whose level spacing fluctuations is expected to be 
described by the gaussian orthogonal ensemble [2§| . 



4- Optical honeycomb potential and graphene 

In graphene sheets, the electrostatic potential that gov- 
erns the dynamics of electrons, the sum of the Coulomb 
potentials of the carbon ions, exhibits the symmetries 
associated to a honeycomb pattern. Of course, in the 
flner details, the optical dipolc potential of ^ and ([5]) 



differs markedly from the graphene potential. In particu- 
lar, the very strong forces that the electrons in graphene 
experience close to the ions have no counterpart in the 
optical lattice, and the interaction between the atoms 
loaded into the optical potential is quite different from 
the electric repulsion between electrons. Nevertheless, 
the common symmetry group implies great similarities 
between the band structures of the two potentials, and in 
the respective parameter regimes where the tight-binding 
approximation is valid, the effective Hamilton operators 
are virtually identical. In particular, experiments made 
with atoms offers new knobs to play with and, with due 
attention to the difference between the two physical sys- 
tems, these observations may deepen our understanding 
about phenomena observed with graphene samples. 

In a very definite sense, the honeycomb potential ^ 
is the simplest of all graphene-type potentials [s^l . Their 
general form is a Fourier sum over the Brillouin vectors, 

V{r) = 5] e^Q ■ with = v*Q . (17) 

QeB 

The various symmetry properties of a honeycomb potcn- 
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tial ensure that the vqs are grouped into sets of closely 
related coefficients. If one coefficient in (|17p is nonzero, a 
whole set of closely related coefficients have correspond- 
ing nonzero values as well. 

Other than the trivial constant solution V{r) = Vq, the 
simplest case is obtained when all coefficients vanish ex- 
cept for the set associated with Wj^^ = Vq and, by conven- 
tion, Vq = 3Vo. This yields the honeycomb potential (O 
with v{r) of dHl). 

III. MASSLESS DIRAC FERMIONS 

A. Band structure in the hopping picture 

In the hopping picture, one envisions the particle as 
hopping from site to site with some quantum mechan- 
ical hopping (or tunneling) amplitude. In the simplest 
situation, all sites have the same energy, only hops be- 
tween nearest-neighbors sites are considered and all hop- 
ping amplitudes take on the same complex value io- The 
one-particle quantum dynamics is then conveniently de- 
scribed using second quantization. In the present situ- 
ation, as we have two different sub-lattices, one has to 
introduce two sets of fermionic annihilation and creation 
operators, one for the A sites, (ajg.,at^), and one for the 

B sites {bj^, bj^), where i and j label the sites in the two- 
dimensional lattices while a stands for the spin index or 
any other pertinent quantum number of the particle. The 
second-quantized Hamilton operator then reads 

(18) 

where means that only nearest-neighbors are in- 

cluded in the sum. The model defined by this Hamilton 
operator accounts for hopping to neighboring lattice sites 
but does not permit a change of the internal quantum 
number a during the hop. We have also included a pos- 
sible energy mismatch e between the A and B sites Q. 
Using the Fourier transform in of the fermionic oper- 
ators, the right-hand side of can be recast into the 
form 

with 

n 

from which wc get the band spectrum 

e±(fc)=±v/e2 + |zj'. (21) 

As expected from the fact that the honeycomb lattice 
consists of two distinct sublattices, we find two bands: 



a conduction band (-I-) and a valence band (— ). These 
bands are here independent of the spin index a, meaning 
that each k G Q accommodates 2<t + 1 internal states 
per subband. Without any real loss of generality, we will 
stick to spin-i fermions in the sequel. As readily checked, 
Zk vanishes when 

^ ^ ^ifc ■ ai ^ ^ifc ■ aa ^ , (22) 

which is solved by the corners K and K' of fl since 
K ■ a2 ^ K' ■ ai = 27r/3. We thus see that the conduc- 
tion and the valence bands are gapped by e, a situation 
typical of a metal when the lattice is filled with particles. 
When there is exactly one particle per site (a situation 
known as half-filling), all levels in the valence band are 
filled at zero temperature, and the Fermi energy Ep (the 
energy of the highest filled level) precisely cuts the en- 
ergy surface at the K and K' points. In this case the 
low-energy excitations of the system can be described by 
linearizing the band spectrum in the neighborhood of K 
and K' . Denoting by q ~ p/h the small displacement 
from either K or A'', the linearization of Zk gives 

|^fe|«^^|q| = ft^'o|ghH^^o, (23) 

where the quantity vq — 3a|to|/(2/i) is called the Fermi 
velocity in the solid-state community. We adopt this ter- 
minology although it is somewhat unfortunate, because 
it has not hing to do with the standard Fermi velocity 
\/2Ep jm^ which depends on the actual mass of the par- 
ticle. 

The dispersion relation now takes on the very sugges- 
tive form 

e±(p) w ±^mlvl+j?-vl (24) 

that is typical of a relativistic dispersion relation with 
particle- hole symmetry. The effective mass m*, defined 
through e = tti^Wq, appears thus as the rest mass of the 
excitations and relates to the energy imbalance of the two 
sub-lattices. The Fermi velocity vq is the analog of the 
velocity of light in relativity. 

The effective Hamilton operator that is derived from 
these considerations and describes the dynamics of the 
excitations around K and K' ^ 

(25) 

where '4>(v) is a 4-componcnt Dirac spinor encapsulating 

- + / 7° \ 

the excitations around K and K' while ?/) = f/jT ( j j 

generates an equation of motion that resembles the Weyl- 
Dirac equation in two dimensions. This is why the name 
Dirac points is given to K and K' (see Refs. 0, @ for 
more details). In this two-dimensional context, the Dirac 
matrices are 7^ = (7^17) = (f^ , icTj; , idj, ) in terms of the 
standard Pauli matrices. 
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FIG. 5: The tight-binding band structure of graphene (in 
units of the tunneling strength \to\) as a function of fe £ SI 
in units of k = \/3fcL- The origin of energy has been chosen 
at the Dirac points and the axis ranges are < 1/2 

and < \/3/3. The bottom contour hnes are hnes of 

constant e / to ■ 



When e vanishes, as is the case in real graphene where 
all lattice sites have the same energy, then e±(fc) = ±|2;fc| 
and the two bands are degenerate at the corners of 
57 where they display circular conical intersections (sec 
Fig. [5]). In the literature, this situation is referred to as 
a semi-mctal or a zero-gap semi-conductor and the cor- 
responding low-energy excitations arc known as massless 
Dirac fermions. The total band width is = 6\to\ and, 
at half-filling, the Fermi energy Ep = 3|<o| (taking the 
energy origin at the lower band minimum) precisely slices 
the energy bands at the Dirac points. Hence the Fermi 
surface reduces to these two points, so that the density 
of states vanishes there see Fig.[H] 



B. Tight-binding approximation 

Mindful of possible experiments, the hopping param- 
eter Iq appears to be an important amplitude to eval- 
uate. We report three different methods for estimating 
its strength |to|- We will start with the familiar tight- 
binding approximation using localised Wannicr functions 
[3lll32j| that are further approximated by Gaussians. We 
will then develop a more accurate semi-classical calcula- 
tion based on an instanton approach (ssj . We will com- 
pare both results to a brute-force exact numerical com- 
putation. 

As a consequence of Bloch's theorem [2^ [s^j , the en- 
ergy spectrum of an atom of mass m moving in the hon- 
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FIG. 6: The noninteracting density of states per unit cell and 
per spin component p{£) as a function of the reduced energy 
£ = _B/|to|. The origin of energy has been chosen at the 
Dirac points. When £ ^ I, then p{£) ~ 2|f |/(\/37r) and the 
density of states vanishes at f = 0, a signature of the semi- 
metal behavior. Note the logarithmic Van Hove singularity 
at \£\ = 1. 



eycomb lattice potential is obtained from 



2m 



Vir) 



(26) 

where we dropped the spin index a which is not essential 
here. The Bloch waves ipnk are given by 



ipnkir) 



Ak ■ r 



Unk{r) 



(27) 



with k € fl, n the band index, and Unkif) is a i3-pcriodic 
function. The latter can be conveniently expanded using 
Wannicr functions [1^, [13, [IHl in accordance with 



(28) 



ReB 



Wannier functions are very useful in describing models 
where particles are localized in space, such as the Hub- 
bard model [36j . They form an orthonormal basis set of 
functions centered at different Bravais lattice sites which 
are copies of the same "seed" functions defined in a given 
primitive cell. The localization properties of the Wannicr 
functions crucially depend on the analyticity properties 
of Unk as a function of k and decay exponentially in the 
simple cases [13, [H, [H, [i^ . 

In the tight-binding approximation, the atoms are as- 
sumed to be sufficiently deeply-localized in the different 
potential wells where they only populate the lowest vi- 
brational levels. Vibrational states in different wells are 
also assumed to have small overlap: the atomic motion 
is thus "frozen" except for the small tunneling amplitude 
between neighboring wells and are then effectively con- 
fined to move in the lowest bands of the lattice. Since 
the Wannicr functions display the same symmetry as the 
local potential structure [4lj, [4^ , the natural idea here is 
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thus to construct tight-binding Wannicr functions from 
hnear combinations of wave functions deeply-locahzed in 
the two potential wells of the primitive cell (the so-called 
atomic orbitals) [13, |4^. This trial wave function ex- 
ploits at best the sub-lattice structure of the honeycomb 
lattice and should give good results at least for the first 
two bands. 

After dropping the band index n, this approach, rem- 
iniscent of the LCAO method (linear combination of 
atomic orbitals) [H, [s^l , leads to the ansatz 



where the quasi-Bloch wavcfunctions 

A 



(29) 



(30) 



essentially live on the type-A sublattice and the type-B 
sublattice, respectively. The sublattice Wannicr func- 
tions Wa(t') and Wnir) are normalized to unity. In the 
present case, we even have w-B{r) = Wj^{—r) due to the 
reflection symmetry of the potential, V{—r) = V{r) [4l| . 
We define the on-site energies as Ea = {walHlwa) {a = 
A, b) and use the parametrization = Eq + A and 
i^B = £^0 ~ A in the following, with Eq the mean on-site 
energy and A half the on-site energy difference. Most im- 
portantly, the sublattice Wannier functions are orthogo- 
nal. However, obtaining their exact expressions is a diffi- 
cult task and one often resorts to simple approximations 
that do not have this property. This is why, in view of 
this very common practical situation, we will consider 
in the following that the Wannicr functions Wa('") and 
Wsir) can overlap. 

Plugging now the ansatz ((29l) - p0)) into ((26l) . and only 
considering coupling between nearest-neighbor lattice 
sites, we get the 2x2 homogeneous linear system 



A^E 



ERk\ 



Zl-ERl -{A + E)J 



0, (31) 



where E = e{k) — Eq and with the matrix entries 

a 

ta = {W^\{n - Eo)\w^J , 



U'a WbJ C 



ifc • Ca 



(32) 



Here Bq = A + Cq is a short-hand notation for the three 
B sites next to the A site. 

Several remarks are in order. First one notes that the 
off-diagonal matrix entries depend on the energy as soon 
as the sublattice Wannier functions overlap. Second, as 
readily checked, the hopping amplitudes ta and E arc 
independent of any energy shift in the Hamilton operator 
and are thus independent of any particular choice for the 



energy origin as one expects. Note also that the values 
of Ej^ and of E^ do not depend on the particular choice 
for point A or point B since Ti. is ^-translation invariant. 
By the same token, the values of ta and of {waIwb^) do 
not depend on the particular choice of A, but B must be 
one of its three nearest neighbors. 

To have a nonzero solution, the 2x2 determinant asso- 
ciated to ((3T|) has to vanish, from which we get the band 
structure. When the overlaps of the sublattice Wannier 
functions is small, {w,^\wb^) <C 1, the band structure is 
very well approximated by 



e±{k)^Eo±^jA^ + \Zk 



(33) 



a form reminiscent of (j21[) . For the honeycomb lattice, 
for which Ti. is yS-pcriodic and invariant under reflections, 
we further have i?A — E^ ~ Eq and A = 0, which implies 
that the effective mass to* of the Dirac fermions is in- 
deed zero. As a consequence, we get the two first bands 
as e±(fc) = Eq ± \Zk\- Furthermore, since V{r) is also 
invariant under 27r/3 rotations about any lattice site A, 
all three tunneling amplitudes ta from A to Bq acquire 
the same value and Z^ of p2|) turns into Zk of pOl) with 



to = y (dr) w:{r){n ~ Eo)wA{r - c) . 



(34) 



where Ti is the differential operator of ([26)) and c is either 
one of the three displacement vectors in p6|) . 



C. Harmonic approximation 

To proceed further one needs an approximation for the 
Wannier functions Wa and Wb- One possibility is to rely 
on the harmonic approximation of the potential wells 
around sites A and B, that is to approximate and Wb 
by the corresponding harmonic ground state wave func- 
tions. We find 



V{ra+r) « -T^o'^V 



""0 „2 



r for a = A, B 



with hwo — S^/VqEji . 



(35) 



where Er = h^k\/ {2m) is the recoil energy of the atom. 
In terms of £ = yjh/ (niLOQ), the familiar length unit of 
the harmonic oscillator, the ground state wave function 
is 



WaIta + r) = WB{rB + r) 



e-hrV^' . (36) 



From this we get Ea = Eb = Eq hujQ and the overlap 
integrals are simply 



\'Wa\Wb, 



exp 




(37) 
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Keeping in mind that Vo ^ Eq ^ in the tight- 
binding regime, {wj^Iuib^) <C 1 and we find from (j34p 



as 




at leading order. However, since the hopping amplitude 
is given by the overlap integral of the localized wave func- 
tions Wa and of two neighboring sites, we see that the 
value of to crucially depends on the tails of these wave 
functions. Wannier functions often decay exponentially 
and, therefore, they cannot be realistically approximated 
by Gaussian wave functions. Hence ([55)) can, at best, 
serve as a rough underestimate [ij. In the next section 
we will derive a reliable and accurate estimate of the tun- 
neling amplitudes in the tight-binding regime by use of 
the instanton method. 



D. Semiclassical estimate 



Using k^^, y/Vo/m, Vq, and ^m/(fc^Vb) as length, 
velocity, energy, and time units, respectively, the 
Schrodinger equation can be conveniently recast into a 
dimensionless form that features an effective Planck's 
constant he (we keep the same symbols for the rescaled 
variables for simplicity), 



v{r)-)jj, 



2Ef 

Vq 



(39) 



with v{r) given by ([5]), here expressed in rescaled units. 
In the tight-binding approximation it is assumed that 
Vq 3> Er, and thus he <^ 1. In this situation, semiclassi- 
cal methods provide very efficient and very accurate ways 
for evaluating dynamical and spectral quantities of inter- 
est. They generally amount to evaluating integrals with 
the aid of semiclassical expressions for the quantum prop- 
agator, derived from its Feynman-path integral formu- 
lation through stationary-phase approximations around 
the classical trajectories [45 1. 

For example, it is well-known that the energy splitting 
between the two lowest energy levels of an atom mov- 
ing in a one-dimensional symmetric double well can be 
accurately calculated using the WKB method [i^ . This 
WKB method can be extended to several dimensions and 
in the sequel we will derive a semiclassical estimate of to 
for the honeycomb lattice using the method proposed in 
[33| . It amounts to evaluating using the classical com- 
plex trajectory (in rescaled units) that connects A and B 
through the classically forbidden region — the so-called 
instanton trajectory. 

Using hioo as an order of magnitude for the vibrational 
level inside a potential well, we see that in the rescaled 
units, this energy is huo/Vo = 3he/V^ ^ 1- So we can 
simply look for the instanton trajectory at zero energy. In 
rescaled units, the hopping amplitude is then expressed 



N 



(40) 



where Sq is the (rescaled) classical action along the zero- 
energy instanton trajectory, and the numerical factor a 
is obtained from integrating out the fluctuations around 
the zero-energy instanton trajectory (see below). 

As the zero-energy instanton fully runs in the classi- 
cally forbidden region, the variables take on complex val- 
ues. For our particular case, the good parameterization 
turns out to keep r real while taking t = it and p = —ip 
purely imaginary with t and p real. Hamilton's classical 
equations of motion in the new variables are just obtained 
from the original ones by flipping V{r) to —V{r). The 
symmetry of the potential dictates that the zero-energy 
instanton trajectory is simply the straight line connect- 
ing site A to B (see Fig. H]) . In the following we calculate 
the instanton between A and A+Ci . Integrating the equa- 
tion of motions, one gets the instanton trajectory in the 
rescaled form to(t) = kLaxo{T)ex with 



tan[7rxo(r)/3] = -VS coth[3V2T/4] 



(41) 



The boundary conditions are xq = 1, io = when 
T —t — oo and xq = 2, xq = Q when r — > oo, meaning that 
the instanton starts at A with zero velocity and ends at B 
with zero velocity, the whole process requiring an inflnite 
amount of time. This is indeed what is expected as both 
endpoints of the instanton are instable in the reversed 
potential picture. Since the energy associated with this 
instanton trajectory is zero, the classical action is simply 

2ka 

So^ ydx|/(x,2/ = 0)|=4V2(^l-^^ « 2.237, 

ka 

(42) 

where f{x,y) is given by ([7]). 

The computation of a proves technically more demand- 
ing. Following [ssj . it is given by the product aia2 with 



ai 



Oi2 



27r y det'[-a2+w2(T) 




det[-a2 



det^a^ 



-UJn 



(43) 



Here w^(t) ~ (5^u)(ro) (a ~ x, y) is the curvature of the 
rescaled potential along the zero-energy instanton tra- 
jectory roir) while ujq is the curvature of the rescaled 
harmonic potential approximation around A; see psp . In 
rescaled units, we have ujq = 3/\/2. The prime in the for- 
mula for ai means that the determinant is calculated by 



excluding the eigenspace of the operator —df 



with 



the smallest eigenvalue. 

The determinants of the differential operators involved 
in the computation of a stem from the linear stability 
analysis of the dynamical flow in the neighborhood of the 
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zero-energy instanton trajectory as encapsulated in the 
monodromy matrix. They can be straightforwardly com- 
puted from solutions of the linear Jacobi-Hill equations 
of degree 2 associated with these differential operators 
[4^. For example, a2 is solved as 



a2 



lim , 

T^oo 



J[T) 



(44) 



where the Jacobi fields J(t) and Jo(t) satisfy the differ- 
ential equations 



i0l{T)J{T) ^ 0, 



dr2 

with initial conditions 



iotMr) = 0, 



(45) 



Jo(-T) = J(-T) = 0, 
io(-T) = j(-T) = l. 



(46) 



The interested reader is referred to [33, Sg] for details. 
We simply give here the final result for the honeycomb 
lattice: 



Oil 



'27^/2 



3.486, a2« 0.449, a « 1.565. 



(47) 

Recasting the semiclassical calculation of the tunneling 
amplitude in units of the recoil energy finally yields 



J-^ w 1.861 — 



3/4 



exp 



-1.582 




(48) 



The same type of scaling laws has been obtained in the 
case of the two-dimensional square optical lattice [4^.[47t. 
In the square-lattice geometry, however, the potential is 
separable and the semiclassical calculation proves much 
simpler as it reduces to using the well-known Mathieu 
equation. 



E. 



Numerical computation of the band structure 



Using Bloch's theorem and plugging ([27|l into ([26]) . we 
get a family of partial differential equations for the MnfcS 
labeled by the Bloch vector fc e fi. After scaling variables 
with the same units as in the previous paragraph, the 
band structure is then extracted by numerically solving 



T-ikUnki.r) = enkUnk{r), 

Hk = -^i-iV + kf + vir) 

for each fc e (expressed now in units of k^). 



(49) 




V3/3 



FIG. 7: Numerically calculated band structure of the two low- 
est energy bands for he = 0.25 at discrete points in the Bril- 
louin zone fi. The same conventions as in Fig. [S]are adopted. 
The value of |to| is determined by requiring that e± = ±3|to| 
at the center F of the Brillouin zone. The similarity with 
Fig. [5] shows that at Vb = 32Er the tight-binding regime has 
already been reached. 



The u„fcS being yB-periodic, they are conveniently 
Fourier expanded in the reciprocal lattice B according 
to 



(r) 



E 



a 



iiQ 



(50) 



The matrix representation of Tik is sparse and banded. 
It is then appropriately truncated and diagonalized such 
that only a small number of coefficients C„q are actually 
significant for the corresponding energy bands. The en- 
ergy bands obtained in this way arc exact and one can 
investigate their dependance on as done in Fig. [7] and 
Fig. [HI 

The essential feature is to realize that the band de- 
generacies at points K and K' are generic and do not 
depend on the actual value of the effective Planck's con- 
stant. Indeed the existence of two degeneracy points in 
the first Brillouin zone for the honeycomb lattice is a gen- 
eral consequence of the lattice symmetries [H, . The 
lattice symmetries are encapsulated in the point group of 
the lattice which is the set of operations that leave fixed 
one particular point of the lattice. The corresponding 
elements are rotations, reflections, inversions and their 
combinations. Combined with ;S-translations, one gets 
the space-group of the lattice. The graphenc point group 
has been analyzed by Lomer [4^ and contains twelve ele- 
ments. In terms of Bloch wave functions ipnk, the lattice 
space-group operations translate into point group opera- 
tions on fc e f2, possibly followed by a reciprocal lattice 
translation to bring back the resulting new wave vector 
in ri. The key point is that degeneracies can only occur 
at Bloch wave vectors which are invariant (up to recipro- 
cal lattice translations) under the action of a nonabelian 
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FIG. 8: Band structure for nearly-free particles moving in a 
weak honeycomb optical potential in units of Vq. The first 2 
levels are plotted as a function of ky/kL at k^/kh = \/3/2, 
so along the vertical edge of Q, from K2 to K'-^, see Fig. [2] 
The solid curves are obtained for he = ^J2Er/Vq = VTO 
and the dashed ones for fte — \/5- As one can see the band 
structure is rather flat in the band centre but the levels curva- 
ture increases when is decreased. The Dirac degeneracies 
in the ground state obtained at the Brillouin zone corners 
are generic and can be inferred from group-theoretic consid- 
erations. Note however that the conical intersections do not 
extend much over the first Brillouin zone when the potential 
is weak but start to spread when is decreased. 



subgroup G of the point group. For the graphene this 
happens at the Dirac points. For example, at corner Ki^ 
beside unity, G is made of two rotations of angles ±27r/3 
about the centre F of the Brillouin zone and three re- 
flexions about the lines connecting F to the three points 
labeled K. This group of order six admits an irreducible 
two-dimensional representation which explains the band 
spectrum degeneracy at the Dirac points. 

This can be nicely illustrated in the weak Vq limit (or 
equivalently when he is large). In this case, the particles 
are quasi-free and the band spectrum can be understood 
in two steps. First, one folds the parabolic dispersion 
relation of the free particle into the first Brillouin zone 
(repeated-zone scheme [l^) and then one couples cross- 
ing levels at Bragg planes by the weak potential. At 
Ki, three plane waves fold with the same kinetic en- 
ergy, namely Ki ~ fci, K2 ~ fc2 and = (see 
Fig. (2) • The weak periodic potential then couples these 
three plane wave states and the coupling matrix elements 
are all identical. The eigenstates of this 3x3 matrix split 
into a singlet and a doublet. When Vq is negative, the sin- 
glet is the ground state which is consistent with the trian- 
gular Bravais lattice obtained in this case ((5 < 0). When 
Vq is positive {S > 0), the doublet becomes the ground 
state and features the tip of the conical intersection be- 
tween the two sub-bands when the quasi-momentum is 
moved away from K, see Fig. [5] 

From the exact numerical calculation, one can extract 




FIG. 9: The hopping parameter \to\ in units of the recoil 
energy Er (crosses) as a function of the inverse of the effec- 
tive Planck's constant he — ^/2Er/Vo as obtained from the 
exact numerical computation. The harmonic approximation 
(dashed curve) and the semiclassical calculation (solid curve) 
of the hopping parameter have been added for comparison 
even if their range of validity is restricted to the tight-binding 
regime he <^ 1. 
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FIG. 10: The hopping parameter \to\ (in units of the recoil 
energy Er) as a function of the inverse effective Planck's 
constant he — ^/TEii/Vo in the tight-binding regime where 
he <ti 1. As one can see, the harmonic approximation (dashed 
curve) is completely off. For example at Vb = 32iJ_R (or 
he = 0.25) |fo| is underestimated by a factor 10 and the dis- 
crepancy gets worse as Vq increases. On the other hand, the 
agreement between the semiclassical calculation (solid curve) 
and the exact numerical computation (crosses) just proves 
excellent. 



the slope of the dispersion relation at the Dirac points 
and then the corresponding tunneling strength |<o| as a 
function of h~-^ , see Fig. [3] Figure [TU] gives the com- 
parison between the exact calculation, the harmonic and 
the semiclassical calculations as a function of in the 
tight-binding regime where ?ie <C 1. As one can see, the 
harmonic approximation is way off whereas the semiclas- 
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e(k) 




FIG. 11: Cut of the linear dispersion approximation along 
Oy at kx = in the first Brillouin zone as compared to the 
actual band spectrum in the tight-binding regime. At half- 
filling, the Fermi energy cuts the band spectrum at the Dirac 
points K and K' . Doping the system away from half-filling 
moves the Fermi energy up or down but the system can still 
be described in terms of massless Dirac fermions provided 
a|q| <^ 2, i.e. provided the change in the Fermi energy is 
much less than the band- width W = 6\to\ itself. By the same 
token, thermal excitations of the system can still be described 
as thermal massless Dirac fermions provided ksT ^ W . 

sical estimate proves excellent. 

F. Reaching the massless Dirac fermions regime 

To access the massless Dirac fermions regime one first 
needs to completely fill the ground state band alone, a 
situation known as half-filling. This is achieved by having 
exactly one fermion with spin state a ~ ±1/2 per Bravais 
cell. Starting from a spin-unpolarizcd cloud of fermions, 
half-filling is thus reached by loading the optical hon- 
eycomb lattice with exactly 2 fermions per Bravais cell, 
corresponding to an average number density p = 1 in the 
tight-binding picture. When this is achieved, the Fermi 
energy slices the band structure at the Dirac points. For 
experiments that study transport phenomena, one would 
also need to subsequently dope the sample away from 
half-filling such that the Fermi energy of the system is 
varied in the linear part of the band structure. 

In a usual experiment, atoms are generally held in an 
external harmonic potential and the optical lattice po- 
tential is superimposed. Reaching half-filling could then 
be done in two steps, first by significantly increasing re- 
pulsive interactions U between fermions through a Fes- 
hbach resonance and then by driving the system into 
the Mott-Hubbard phase with one fermion per site as 
done in [sO, HH. Then setting U to zero again should 
maintain the system at number density p = 1. Obvi- 
ous candidates for such experiments are Potassium-40 
as well as Lithium-6 atoms [13, HH, [13 • the exter- 
nal trap, the Mott insulator appears first where the lo- 
cal filling is approximately one atom per site and one 



needs to ensure that adding more atoms (or increasing 
the chemical potential /i) does not favor the appearance 
of the doubly-occupied Mott phase. This will be the 
case for very strong repulsion U ^ /i,to,fcBr in which 
case one expects the entire centre of the trap to con- 
tain a Mott insulating phase with single occupancy and 
negligible thermally-activated doubly-occupied sites. In 
the case of the honeycomb lattice, starting from a spin- 
unpolarized sample, it is known that half-filling is reached 
for Uc ~ 5to, the atoms displaying at the same time an 
anti-ferromagnetic order [l^. Note that Uc ~ W, where 
W = 2Ef = 6\to\ is the band- width. 

Doping the system could be done in the following 
way. The external harmonic confinement (with angu- 
lar trap frequency f2t) defines a characteristic length 

C = ^2|io|/(TOrit) over which the energy is shifted by 

precisely the tunneling energy to [H, [13, [H| • This length 
defines the distance over which one given lattice site is 
coupled by tunneling to its neighbors. In turn, having 
loaded Np fermions into the trap, one can define a char- 
acteristic filling factor through p — Np{a/Qi^. Varying 
I to I by changing the lattice potential height Vb or tighten- 
ing/loosening the trap by changing l^t would thus allow 
to tune /5 in a controlled way and hence to dope the sys- 
tem. 

For the conical intersection at the Dirac points to sig- 
nificantly spread over the Brillouin zone f2, one needs to 
reach the tight-binding regime where Vq is large enough 
(typically Vq > 10£^h will do). Inspection of the Tay- 
lor expansion of ([^0]) then shows that it is sufficient to 
have |q|a ^ 2 (q being the small displacement from a 
Dirac point) for the band structure to be well approxi- 
mated by a linear dispersion relation around the Dirac 
points. The available energy range Ai? is thus set by 
the band- width itself, namely Ai? <^W. So tuning the 
filling factor away from half-filling and residual thermal 
fluctuations will keep the system in the massless Dirac 
fermions regime provided /i, fc^T <C W (Fig. [TT|) . For 
example, at Vo = S2Eii, the temperature constraint, as 
derived from (gH]), is T < Tr/ 50 whereas it is T < Tn/2 
at Vq = lOEji. There is thus room left for reaching the 
massless Dirac fermions regime within the current state- 
of-art cooling technology. 



IV. ROBUSTNESS OF THE MASSLESS DIRAC 
FERMIONS 

As the very existence of the massless Dirac fermions 
regime rests on the two conical degeneracies in the band 
structure, one may wonder if this regime would resist 
imperfections of the system. Indeed the argument we 
gave to explain the conical degeneracies relied on group- 
theoretic arguments which were specific to the hexagonal 
symmetry of the honeycomb structure. In practice, it is 
impossible to control the laser configuration to the point 
where all intensities and alignment angles would all be 
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exactly equal. Such imperfections in the system would 
obviously break the hexagonal symmetry and one could 
think that the Dirac fcrmions would just be destroyed. 
In fact, as we will see shortly, massless Dirac fermions 
are quite robust and survive small imperfections that are 
easily within experimental reach. 



A. Imbalanced hopping amplitudes 

To understand why massless Dirac fermions are robust, 
we will start by analyzing the case of imbalanced hop- 
ping amplitudes as done in (SGj . For real graphene, this 
would correspond to stretching the graphene sheet. In 
this case, the tight-binding band structure is given by 
e±{k) = ±\Zk\, where Zk is defined in ^0^. The de- 
generacies are found at points ko & canceling ^ 0. 
This condition boils down to sum up three vectors to zero 
in the two-dimensional plane, see Fig. [T^l As such, a so- 
lution is only possible provided the hopping amplitudes 
satisfy one of the norm inequalities given by 



<2 - \t.\ <\h\< u 



t3 



(51) 



and cyclic permutations. If this is the case, defining the 
angles ipi,2 = argt2,3 — Stvgti, the Dirac points solve 



cos(fc£) -ai — ipi) 



cos(fcD-a2 — 1^2) 



2\tlt2\ 





1 


4' 




2\ht3 





(52) 



subject to the condition 

I sin(fci3-ai - (fx) + \h\ sin(fc£)-a2 - ^2) = 0. (53) 



t2 



We find the important result that the system self-adapts 
to changes in the hopping amplitudes by shifting the 
Dirac points away from the corners of the Brillouin zone 
until the norm inequalities (|5ip break and degeneracies 
disappear. Thus, provided the hopping imbalance is not 
too strong, the massless Dirac fcrmions do survive im- 
perfections in the system and the hexagonal symmetry 
breaking. 

We illustrate this important feature in the simple case 
of only one imbalanced hopping amplitude, namely t\ = 
ito, t2 = ^3 = to. We further choose 7 real and < 
I7I < 2 for the Dirac points to exist. We then find two 
Dirac points D-y and D'^ given by ko = —k'jj = ipo{b2 — 
bi) where (po G [0,1/2] solves cos(27riy9o) = ^7/2. This 
means that the two Dirac points and D'y move along 
opposite paths in the Brillouin zone fi. The fact that 
Dirac points always come in by pairs of opposite location 
in Q is generic [57[- When 7 is increased from to 2, 
starts at fco = (3^^/4)6^ for 7 = 0, then moves along 
axis Oy and reach corner Ki at 7 = 1. Note that when 
7 — > 0, the physical situation is that of weakly coupled 
"zig-zag" linear chains. For 7 > 1, leaves fl but 




FIG. 12: The condition Zk = is equivalent to cancel the 
resultant vector u of three vectors, each with length \t„\ and 
polar angle a„ = fc ■ c„ -f arg(t„). There will always be a 
solution provided one of the norm inequalities ||t2| — |i3|| < 
1^1 1 < 1^2! + 1^3 1 (a-nd cyclic permutations) is satisfied. 



a translation in reciprocal lattice brings it back on the 
vertical edges of fi (technically we get two copies of the 
same point). reaches the middle of the vertical edge 
at 7 = 2 where it merges with D'^ into a single Dirac 
point, see Fig. [131 Interesting physics occurs at 7 = 2 
in connection with the quantum Hall effect [l3, Ell . As 
soon as 7 > 2, the degeneracy is lifted and the massless 
Dirac fermions do not exist anymore. For negative 7, D-y 
and Di^ move back from ±(3/sL/4)ey to the centre F of 
the Brillouin zone where they merge and disappear, see 
Fig. [131 The fact that Dirac points can only merge at 
the centre and mid-edge points of fl is also generic \5l\ . 

Hence, far from being a nuisance, we see that con- 
trolling the hopping amplitude imbalance proves an in- 
teresting way of exploring the massless Dirac fermions 
physics under different circumstances by moving around 
the Dirac points in the Brillouin zone. 



B. Optical lattice distortions 

The previous discussion concentrated on the impact 
of imbalanced hopping amplitudes irrespective of the 
change of symmetry of the lattice potential. We will 
now analyze these lattice distortions in more detail and 
give quantitative estimates about the experimental de- 
gree of control which is required to target the massless 
Dirac fermions regime. We will consider in-plane laser 
beams with different (positive) strengths En ~ s„i?o and 
with respective angles away from 27r/3, see Fig. [TH It is 
important to note that we will always stick to imperfec- 
tions which are compatible with a two-point Bravais cell. 
They will only induce distortions of the hexagonal spa- 
tial structure of field minima but without breaking this 
pattern. 

The new optical lattice potential is now given by 
V'{r) = Vo \f'{r)\ with the new total dimensionless field 
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D'^ = D2 




■.I D-2 = D' 



FIG. 13: When the three hopping amphtudes t„ are unbal- 
anced, the Dirac points are shifted in the Brillouin zone Q 



and disappear when the norm inequality ^2 



m < h < 



\t2\ + \t3\ is no longer satisfied. We depict here how the Dirac 
points and Diy move in when only one hopping ampli- 
tude is imbalanced, namely ti ='yto and t2=t3=to. Points D-y 
(thick path) and D'^ (thin path) move along opposite paths. 
Increasing 7 from 0, point Dj starts at Do and moves upward. 
It reaches point Ki at 7 = 1 (balanced amplitudes case) then 
moves along the vertical edge of n where it reaches its middle 
point D2 at 7 = 2. The Dirac points cease to exist when 
7 > 2. For negative 7, moves downward from Dq (dotted 
thick path), reaches the zone center F for 7 — —2 and then 
ceases to exist for 7 < —2. 



amplitude 

/'(r) = si + 8-2 cxp(-i6' • r) + S3 cxp(ib2 • r). (54) 

Here the b'^ (n=1.2) feature the new reciprocal lattice 
basis vectors. They define in turn a new set of Bravais 
lattice basis vectors giving rise to a new primitive 
diamond-shaped cell S'. Unless the angle mismatches 
vanish, the new Bravais and reciprocal lattices are no 
longer hexagonal but oblique with no special symmetry 
except for inversion. As a consequence, the new first 
Brillouin zone Q' is still a hexagon but no longer a regular 
one. 

Since we assume a two-point primitive cell, the minima 
of the new optical potential still identify with zeros of 
/'(r). Similarly with the case of imbalanced hopping 
amplitudes, we find two solutions if and only if the field 
strengths s„ satisfy one of the norm inequalities |s2 — 
S3I < si < S2 -f S3 (and cyclic ones). In this case the 
minima are given by 



cos(6' • r) 
cos(62 • r) 



2siS2 
„2 



2siS3 



(55) 



subject to the condition S2 sin(6' • t) = S3 sin(62 • f )■ 

In the following we will examine separately the effect 
of strength imbalance and angle mismatch. 




FIG. 14: (a) [Color online] (a) The asymmetric in-plane 
3-beam configuration. Three monochromatic and linearly- 
polarized laser beams with wave vectors k„ interfere with dif- 
ferent strengths E„ = SnEo {n = 1, 2,3). The respective an- 
gles depart from 27r/3. (b) Distorted optical lattice obtained 
with = = 5 X 10"^ and si = 1, S2 = 1.03, S3 = 0.97. For 
weak enough distortions, the primitive diamond-shape cell E 
still contains two field minima as evidenced in the plot. 



1. Critical field strength imbalance 

To give an estimate of the critical field strength im- 
balance beyond which the Dirac points cannot survive, 
we consider the simple case of only one imbalanced laser 
beam and no angle mismatch, namely O2 = 0^ = 0, 
si = 1 + r/ and S2 = S3 = 1. In this case the Bravais 
lattice, the reciprocal lattice, the primitive cell S and 
the Brillouin zone are not modified. The new optical 
potential V'{r) ~ Vov'{r) reads 



v'{r) = v{r) -f 2ri5v{r) + ri{ri H 
Sv{r) = cos(6i-r) + cos(62-'"), 



2), 



(56) 



where v{r) is given by ([5]). Note that when only one field 
strength is imbalanced, the corresponding potential still 
displays a reflection symmetry. In the present case, it 
is the Ox- reflection symmetry because V'{r) is invariant 
under the exchange bi <-^- 62- Requiring now that the 
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FIG. 15: [Color online] Slightly distorted lattice obtained with 
vanishing mismatch angles and one imbalanced field strength, 
namely si = 10/9 and S2 = S3 — 1. In this particular case the 
hexagon of field minima is slightly squeezed along the horizon- 
tal axis Ox and the vectors c'„ connecting a given minimum 
to its three nearest-neghbors have now diff'erent lengths. In 
the situation depicted |c2| = |c3| 7^ |ci |. In turn, due to the 
reflection symmetry about Ox, the tight-binding hopping am- 
plitudes satisfy |i2| = |i3| 7^ 



primitive cell E exhibits two field minima imposes — 1 < 
77 < 1. Their positions in S are given by r' ^ = (^a,b (0.1 + 
a-z) with cos(27r(pA,B) = ^ (1 + ??)/2. Their mid-point 
r' = (r^ -|- T"g)/2 = (ai -I- a2)/2 is a saddle point and 
defines the potential barrier height to cross to go from 
A and B in E. One finds V/ = (77 - l)^Vb. 

As a whole the field minima organize in a hexagon 
which is stretched (77 negative) or compressed {rj positive) 
along Ox, see Fig. (TSl As a consequence two of the three 
new vectors joining one minimum to its three nearest- 
neighbors will have equal length. In the present situation 
we get IcjI = \c'^\ |c'|. The potential barrier height 
Vg' to cross to go from A to B along c'2 and Cg is given by 
the corresponding saddle points located at the middle of 
the edges of E. One finds V^' = (77 + I^Vq. 

Now, when 77 is increased from 0, the minima move 
closer along c[ and move away along c'2 and Cg. At the 
same time, the potential barrier along c'l is lowered 
and the the potential barrier V^' along and Cg is in- 
creased. As a net effect, in the tight-binding picture, we 
expect the tunneling amplitude to increase while 1 
and \t3 \ decrease. We get the opposite conclusion when 77 
is lowered from 0. Since the potential is invariant through 
bi <-!■ ^2, we further have \t2\ = Its] and we recover the 
case of one imbalanced hopping amplitude analyzed in 
the previous section. 

One could try to derive a semiclassical expression of 
the tn as a function of 77 using the instanton method but, 
actually, such a tedious calculation proves unnecessary, 
at least when 77 is small. Indeed, by inspection of the 
semiclassical expression (|40|). we expect the ratio |ti/i2| 
to scale as exp(A5(7;)/?ie) at leading order, where AS'(77) 



is the action difference between the two instanton trajec- 
tories linking sites A and B along and c'l respectively. 
For small enough 77 we expect AS'(77) to grow linearly with 
77, the slope being positive since the ratio |ii/i2| should 
increase with 77. The Dirac degeneracies disappear when 
this ratio is 2 (see previous section), thus we get the semi- 
classical prediction that this will happen when 77 cx fte- 
This result can also be inferred by saying that the Dirac 
points will disappear as soon as the perturbing potential 
277i5I^(r), see (|56| . strongly mixes the unperturbed states. 
This will happen when the corresponding coupling energy 
equals the mean level spacing of the unperturbed system 
and we get back to the prediction 77 cx fi,e- 

To check our semiclassical prediction we have com- 
puted, for each value of the effective Planck's constant 

he, the ground state and first excited-state levels for dif- 
ferent values of 77 and we have extracted the correspond- 
ing critical value 77c for which the Dirac degeneracies are 
lifted. Figure [TBI gives an example of the band structure 
obtained at = 1/V40 « 0.16 for 77 ranging from to 
0.054. We have then plotted 77c as a function of h^, see 
Fig. [T71 We have fitted the data with the quadratic fit 
function ahe + f3hl and found a w 0.1074 and (3 0.0624 
enforcing the very good agreement obtained with our 
linear prediction in the semiclassical regime ?ie <C 1. 
The quadratic correction could certainly be inferred from 
semiclassical higher-order corrections. 

We would like to emphasize at this point that increas- 
ing or decreasing 77 from is not symmetrical. When 
77 is decreased from 0, the Dirac degeneracies are pre- 
dicted to disappear when |ii/i2| 0. However the best 
that we can do is to let 77 —1. This unfortunately 
means that one laser beam is almost extinguished and 
the situation is more that of very weakly coupled one- 
dimensional chains, a situation we postpone to future 
study as it proves interesting for high-Tc superconduc- 
tivity [5^ . We thus see that decreasing slightly 77 from 
does not harm the Dirac degeneracies. They move inside 

but do survive. By contrast, increasing slightly 77 from 
docs destroy the Dirac degeneracies as soon &s rj ^ hp. 

As one can see from the plots, the tolerance about 
the intensity mismatch of the laser beams increases with 

hf, , or equivalcntly when the optical lattice depth Vq de- 
creases. On the other hand, as we already saw, the Dirac 
cones do not extend much over the Brillouin zone if Vq is 
too small. So there is a trade-off to make. The situation 
is however really favorable since the intensity mismatch 
tolerance is already in the 10% range for Vq ^ lOiJjj. 
This means that the massless Dirac fermions prove quite 
robust and should be easily accessed experimentally. 



2. Critical in-plane angle mismatch 

We now estimate the critical angle mismatch when all 
laser beams have the same intensities (si — S2 = S3 = 1). 
We see from ([5^ that the new optical potential still dis- 
plays the exchange symmetry b'^ and thus a re- 
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0.0015 r. 




FIG. 16: The band diagram for the two lowest levels as a 
function of rj for Vb = 80Er {He ~ 0.158). The bands are 
plotted along the vertical straight line joining the Dirac points 
K2 and K-j of the balanced situation, see Fig. (2] The origin 
of energy is fixed at the Fermi energy for a half-filled band 
and all bands have been shifted such that the upper and lower 
bands intersect at zero energy difference. 
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FIG. 17: The critical laser strength imbalance rjc at which 
the Dirac degeneracies ar e lifted a s a function of the effective 
Planck's constant fte = ^J2Er/Vo. The solid line corresponds 
to a quadratic fit of the numerical data. The linear coefficient 
is Q « 0.1074 while the quadratic one is /3 « 0.0624. As one 
can see our numerical results are in good agreement with our 
semiclassical prediction rjc oc he. The degree of control of the 
intensity imbalance of the laser fields gets more stringent as 
the optical lattice depth Vo is increased. Nevertheless, at al- 
ready Vo = 2QEr {he ~ 0.3), the laser intensities should all 
equal within 8% which does not sound particularly demand- 
ing. 



flection invariancc with respect to their bisectrix. In the 
following we stick to the simple case where 63 = —62 = 9 
and is small. In this case both the Bravais lattice, the 
reciprocal lattice, the Brillouin zone O and the diamond- 
shaped primitive cell S get modified. The new reciprocal 
basis vectors turn out to be 6' = bi + ^61, = ''2 + ^^2 
where 8bi = (6»/V3) ba and ^62 = (6I/V3)bi. Since the 
exchange symmetry bi <-> ba is again preserved, the new 
potential continues to display the Ox-reflection invari- 
ancc. Figure [TH] gives a plot of the new potential struc- 




FIG. 18: [Color online] Distorted lattice obtained with bal- 
anced field strengths Sn = 1 and angle mismatch 63 = —62 = 
— tt/IO. In this particular case the hexagon of field min- 
ima is stretched along the horizontal axis Ox and the vec- 
tors connecting a given minimum to its three nearest- 
neghbors have now different lengths. In the situation depicted 
|c2| = |c3| 7^ |c'|. In turn, due to the reflection symme- 
try about Ox, the tight-binding hopping amplitudes satisfy 

|t2| = |f3|/|il|. 



ture for 6 = —tt/IO. 

This situation boils down again to the case of one im- 
balanced tunneling amplitude. Indeed, the angle between 
the b[ and b'2 decreases when 9 is increased from 0. In 
turn the angle between the corresponding increases 
and the hexagon structure made by the A and B minima 
get compressed along Ox. The opposite conclusion holds 
when 6 is decreased from 0. We get again the situation 
where [is] = l^al \ti\ and \ti/t2\ > 1 when 9 > 
and vice-versa. Like for the field strength imbalance, the 
situations 9 > and 9 < are not symmetric. The mas- 
less Dirac fermions prove more sensitive to closing the 
angle between the b^, so for 6*3 = —92 = 9 > because 
1^1/^2! then increases and the threshold |ti/i2| = 2 is 
more rapidly hit. This is the situation we explore. 

Applying the same reasoning as before, we thus predict 
the critical angle mismatch beyond which the masslcss 
Dirac fermions are destroyed to scale as 9c oc he. Again, 
to get 9c as a function of fig, we numerically compute the 
band structure at a given he for different in-plane mis- 
match angles 9 and then extract the value 9c for which 
the Dirac degeneracy is lifted. We then repeat the pro- 
cedure for different he. As one can sec, our prediction 
is in very good agreement with the numerical calcula- 
tions, see Fig. 1191 and well supported by a quadratic 
fit. As 9c increases with he, there is a trade-off to make 
between reaching the tight-binding regime where Vq is 
large and achieving an experimentally reasonable angle 
mismatch tolerance which requires Vq to be small. The 
trade-off turns out to be a favorable one since already for 
Vo = 2QEjj {he ~ 0.3), one gets a tolerance of about 5° 
on the laser beams alignment. We expect the same type 
of scaling for small out-of-planc angle mismatches. Fur- 
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FIG. 19: The critical angle mismatch 9c (in units of tt) beyond 
which the Dirac degeneracies disappear as a function of the 
effective Planck's constant he = y^2ER/Vo. The dashed line 
corresponds to a quadratic fit of the numerical data. The lin- 
ear coefficient is 0.109 while the quadratic one is —0.0577. As 
one can see our numerical results are in good agreement with 
our semiclassical prediction 9c oc he. The degree of control of 
the angle mismatch gets more stringent as the optical lattice 
depth Vo is increased. Nevertheless, at already Vb = 20Er 
{he « 0.3), the angle mismatch should be less than 5° which 
is not particularly demanding. 



thermore, when several small distortions combine, their 
effects should add up and thus the critical imperfection 
threshold should still scale with he- 

As an overall conclusion we see that massless Dirac 
fermions are quite robust to moderate lattice distortions. 
Demonstrating them in an experiment should not be par- 
ticularly demanding in terms of the control of the laser 
configuration. 



imposing three independent standing waves, of the same 
wavelength and with equal intensity, whose wave vectors 
form the trine of Fig. [T] As a consequence of the 

incoherent superposition, the t replacement of ^ is not 
available, and the r replacement alone cannot remove all 
three phases of the standing waves. One can, however, 
shift r such that the three phases are the same, and then 
one has an intensity pattern proportional to the right- 
hand side of ((57)) . 

Most of the hexagon structure of Fig. 5] remains un- 
changed by this modification: lattice sites A, B, C con- 
tinue to be the locations of local minima and maxima, 
whereas the saddle points S acquire new positions on the 
. . . ABCABC. . . lines. 

Figure [^ confirms that, for small ip values, the minima 
of the honeycomb dipole potential are still organized in 
a hexagonal pattern but we now have different potential 
depths at sites A and B. The potential energy mismatch 
is 2e w 8Vb|(/3|/V3. In view of (HIl) and this means 
that the Dirac fermions acquire a mass to* oc or, in 
other words, that the Dirac degeneracies are lifted. The 
possibility of fine-tuning the mass of the Dirac fermions 
through the parameter ip is an interesting experimental 
knob to play with. 

Increasing |(^| further, one can also see that, for the 
particular values |(^| = tt/6 and 7r/2, the three sublattices 
of saddle points merge into a single triangular lattice, 
which coincides with the A, B, or C lattice, respectively; 
see Fig. 1201 This merging of a potential minimum or 
maximum with three saddle points, leads to a peculiar 
third-order saddle point. For (p — 7r/6, say, the S sites 
merge with the B sites and we have 



^ cos (6a 



tp—Tr/6 



(58) 



C. Inequivalent potential wells 

We finally briefly mention how to distort the optical 
lattice in a systematic manner as it allows for an experi- 
mental control of the mass of the Dirac fermions as well 
as for a continuous switch from a honeycomb lattice to a 
triangular one. 

In Sec. Ill B 41 we observed that the honeycomb poten- 
tial ([6]) is the simplest of all graphene-type potentials, 
characterized by choosing vq and Vb^ real (in fact, pos- 
itive) while putting all unrelated coefficients in p7|) to 
zero. Now, letting Vbi to acquire a phase (p, such that 
e~'^^Vbi is positive, will break the reflection symmetry of 
the honeycomb potential [s^l . 

In the r-dependent part of the dimensionless potential 
([8]), this phase ip is introduced by the replacement 

3 3 

^ cos{ba • f ) ^ ^ cos{ba ■ v + ip) , (57) 

a—1 a—1 

where 63 — —61 — 62- This can be implemented by super- 



for |r — ^ K ^, hence a cubic saddle point rather 
than the usual quadratic saddle point. 

An unpolarized ultracold gas of spin-i fermions loaded 
into such a potential at half-filling would lead to two 
fermions per well. By driving the system through attrac- 
tive interactions, one could even get a Mott insulator of 
fermion pairs. By switching off all interactions and set- 
ting (fi = 0, one should be able to study oscillations of 
atoms between the A and B sublattices. We will analyze 
this situation in a follow-up paper. 



V. CONCLUSION 

Motivated by the vivid field of graphene physics, we 
have explained and analyzed how to reproduce massless 
Dirac fermions by loading ultracold fermions in an opti- 
cal lattice with honeycomb structure. We have described 
the two-dimensional laser configuration that gives rise 
to an optical potential where field minima are organized 
in a honeycomb structure (with lattice constant a) and 
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FIG. 20: For various values of the phase parameter of (|57[) . 
the plot shows the potential energy along a . . . ABCABC. . . line 
in Fig.O The top plot, for ip = 0, repeats the bottom-left plot 
of Fig. |4]for reference. The degeneracy between sites A and B 
is lifted for the small ip value of ip — 7r/24, the saddle points 
have moved closer to the B sites, where we continue to have 
local minima. In this situation the Dirac fermions acquire a 
mass m, oc \ip\. When = 7r/6, the saddle points S coincide 
with the B sites, and we have cubic saddle points there. Fi- 
nally, in the bottom plot, we have <p = n/S and the saddle 
points are halfway between adjacent B and C sites, with po- 
tential maxima at both of them. Except for a displacement, 
the potential in the bottom plot is the negative of the poten- 
tial in the top plot, and thus identical with the honeycomb 
potential (|6]) for red rather than blue detuning of the three 
running wave lasers. For ease of comparison, the potential 
constants are adjusted such that the maxima and minima are 
ai V = and V = 9Vc), respectively, for all Lp values. 

we have thoroughly detailed the corresponding crystallo- 
graphic features. The behavior of atoms propagating in 
such an optical potential in the tight-binding regime is 
in one-to-one correspondence with the behavior of elec- 
trons propagating in a graphene sheet. The ground state 
and first-excited levels of the band structure exhibit two 
conical degeneracies located at the corners of the first 



Brillouin zone, as dictated by symmetry arguments. In 
the neighborhood of these degeneracies, the band spec- 
trum is linear. 

When the lattice is loaded with fermions at half-filling, 
the Fermi energy slices the band structure at these de- 
generacy points, known as the Dirac points. Around 
half-filling, the tight-binding Hamilton operator can then 
be recast in a form reminiscent of the relativistic Weyl- 
Dirac Hamilton operator and featuring so-called massless 
Dirac fermions. The important parameter driving the 
dynamics turns out to be the hopping amplitude to be- 
tween nearest-neighbors sites as it gives the band width 
W = 6\to\ and the "Fermi velocity" vq = 3a\to\ / {2h) . We 
have derived a semiclassical expression for \to\ in terms 
of the e ffective P lanck's constant of the problem, namely 
he = y/2Eii/Vo (with Vb the optical potential strength 
and En the recoil energy) and have compared it to an 
exact numerical calculation of the band spectrum. From 
this we have derived quantitative experimental criteria 
(such as the required initial temperature of the atomic 
gas) to reach the massless Dirac fermions regime. 

We have also examined the robustness of the massless 
Dirac fermions to imperfections of the laser configura- 
tion (field strengths imbalances and angle mismatches). 
Massless Dirac fermions turn out to be quite robust as 
the equality of the beam intensities should be controlled 
within the few percent range while the respective beam 
angles should equal 2tt/3 within the few degrees range. 
By appropriately controlling these lattice distortions, one 
can even control and move the Dirac points in the Bril- 
louin zone. Lastly, we introduce an irremovable phase 
to the honeycomb potential, hence lifting the degeneracy 
between two sublattices. In turn, the Dirac fermions ac- 
quire a mass proportional to this phase. We also briefly 
mention the peculiar properties of saddle points and the 
possibility to study oscillations of atoms between two 
sublattices as a consequence of this irremovable phase. 

As an overall conclusion, mimicking graphene physics 
with ultracold fermions is within experimental reach. For 
non-interacting fermions, one could think of implement- 
ing transport experiments (in the presence of disorder or 
not). For example, by rotating the whole honeycomb lat- 
tice around a given axis perpendicular to the latti ce p lane 
or by implementing the scheme proposed in |62{ . one 
would mimic effective magnetic fields able to reproduce 
the quantum Hall effect situation. In the rest frame of 
the atoms, the centrifugal effects arc described by a fic- 
titious vector potential. This leads to Landau levels and 
paves the way to physical effects analogous to the quan- 
tum Hall effect. The possibilty to move the Dirac points 
in the Brillouin zone even offers new physical effects to 
test 

Interacting systems on a lattice prove also particu- 
larly interesting as they can be mapped (at least for 
strong interactions) on Heisenberg models and thus of- 
fer ways of exploring quantum magnetism (63j . In the 
case of the honeycomb lattice, quantum phase tran- 
sitions are predicted to occur when the interaction 
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strength \U\ is strong enough. For repulsive interac- 
tions, quantum Monte-Carlo calculations pred ict anti- 
ferromagnetic order to occur at half-filling [13 • For at- 
tractive interactions, mean-field calculations have started 
to analyze the BEC-BCS crossover and predict a semi- 
metal/superconductor transition [l^. Recent Monte- 
carlo studies have even started to analyze this BEC-BCS 
crossover [64| and one can expect an increase of such 
studies in the near future. Very recently, implementa- 
tions of massless Dirac fermions in square lattices have 
been proposed [1^, [6^. The situation seems thus ma- 
ture for an experimental effort towards loading ultracold 
fermions in a honeycomb optical lattice. 
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